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Abstract. We construct explicit invariant measures for a family of infinite 
products of random, independent, identically-distributed elements of SL(2, C). 
The matrices in the product are such that one entry is gamma-distributed along 
a ray in the complex plane. When the ray is the positive real axis, the products 
are those associated with a continued fraction studied by Letac & Seshadri [Z. 
Wahr. Verw. Get. 62 (1983) 485-489], who showed that the distribution 
of the continued fraction is a generalised inverse Gaussian. We extend this 
result by finding the distribution for an arbitrary ray in the complex right- 
half plane, and thus compute the corresponding Lyapunov exponent explicitly. 
When the ray lies on the imaginary axis, the matrices in the infinite product 
coincide with the transfer matrices associated with a one-dimensional discrete 
Schrodinger operator with a random, gamma-distributed potential. Hence, 
the explicit knowledge of the Lyapunov exponent may be used to estimate the 
(exponential) rate of localisation of the eigenstates. 



1. INTRODUCTION 

Let {-ArajneN denote a sequence of independent random matrices identically dis- 
tributed according to a probability measure /i on GL(2,R). The problem of deter- 
mining the asymptotic behaviour of the product 

(1.1) Mn ■= A1A2 ■ ■ ■ An as n ^ CX), 

plays a crucial role in the theory of products of random matrices and its applications, 
especially in mathematical physics [HlC]. 

The rate of growth of this product can be quantified by its Lyapunov exponent 

(1.2) A := lim -E{\n\Mn\), 

n— >oo n 

where | • | denotes some matrix norm. 

The Lyapunov exponent exists whenever E(log^ l^il) < co- The Furstenberg- 
Kesten theorem (Jl], [HI P- H) generalises the classical strong law of large numbers 
to the case of non-commuting random products and states that 

A= lim — ln|A^„| almost surely. 

n— >oo n 

The case of unimodular matrices (i.e. An G SL(2,R)) is of particular interest. In 
this case, under the natural additional assumption of noncompactness^ and strong 
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irreducibility,'^ Furstenberg's theorem (^Sl) cf- |E] P- 30), asserts that the Lya- 
punov exponent is strictly positive. Moreover, there exists a unique, continuous, 
//-invariant measure v on the projective hne P (M^) — that is, a measure that is 
invariant under the projective action of matrices drawn from the distribution fj, (cf. 
p. 30). 

The calculation of the Lyapunov exponent involves this //-invariant measure, but 
there are remarkably few non-trivial cases where v has been found explicitly; three 
well-known examples will be discussed presently (others are given in j^) . The 
prominent feature shared by these examples is that the invariant measure is found 
by considering a random continued fraction derived from the projective action of 
the relevant matrix ensemble. 

The first example dates back to Dyson's work on the disordered chain problem. 

Example 1 (Dyson [El)- The disordered chain is modeled by a system of harmonic 
oscillators coupled by linear forces. A physical realisation is obtained by considering 
a sequence of N particles joined by elastic springs obeying Hooke's law. Denote 
the mass and the displacement of the nth particle from its equilibrium position by 
m„ and Xn respectively, and let kn be the elastic modulus of the spring between 
the nth and {n + l)th particle. Then the equation of motion takes the form 

TO„ Xn = kn{Xn+i - Xn) + fc„_i(a;„_i - X„), 1 < n < iV . 

By introducing additional variables, it is straightforward to express this as the 
first-order system 



where 



= ^/a^Un+l - y'ttn-l Un-1, 1 < n < 2A^ - 1, 



kn kn 
a2n-l , a2n '■— 

TO„ rrin+i 



In matrix notation, 

u = Ju, 

where is the tridiagonal (Jacobi) matrix 

/ 

-^al 

-V^ ^ 



J = 





V 



'y/a2N-3 









\ 





^Ja,2N-2 

/ 



Dyson studied the spectral problem for J when the a„ are independent, identically 
distributed random variables. It turns out that the spectral properties of J (e.g. 
the eigenvalue density function) can be deduced from the so-called characteristic 
function of the chain 



ln(l + x)vt[d,x\ 



^There is no finite union of proper subspaces V with A-\V = V, for all reaUzations of 
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where vt is the distribution of the random continued fraction 

ait 



(1.3) 



1 + - 



1+ 



1 + 

To iUustratc his approach, Dyson elaborated the particular case where the a„ are 
gamma-distributed, i.e. for every Lebesgue-measurable subset S of R+, 



Pr(a„eS')= / 7p,s(a)da, 
Js 

where 

(1.4) 7p,,(a) -^^flf-^e^^ , a > 0; p, s > . 

Then the probability density function of the random continued fraction H1.3|l is 
given explicitly by 



-*(dx)=C,,,-^exp(-^) 



dx, 



{l + xY 

where Cp^s,t is a normalisation constant. 

In terms of products of random matrices, Dyson's continued fraction l|1.3|l cor- 
responds to the case where 

. _ (O a^t^ 

in the product (|1.1|) . The distribution Vt of the continued fraction is the /^-invariant 
distribution associated with this product. 

One of the present paper's contribution is the calculation of the distribution 
of a random continued fraction that arises in another important physical model, 
namely the discrete Schrodinger equation with a random potential. In this regard, 
Lloyd's model (cf. [201) ED) with a Cauchy-distributed potential in one spatial 
dimension is possibly the best-known example where the invariant measure and the 
corresponding Lyapunov exponent have been found explicitly. 

Example 2 (Lloyd's model 6 , p. 35). Let z = a; + ij/ e C with t/ > and denote 
by Cz the distribution of the random variable x + yc, where c is Cauchy-distributed, 

Pr(ce5) = i/-^ 

for every Lebesgue-measurable subset S of M. Set 

An-.^i^'^ "p^j eSL(2,R), 

where the a„ are independent and Cz-distributed. The invariant measure v is then 
given by the distribution C„, where u and z are related by z = u -I- l/u. It follows 
easily that 

A = In \u\ . 

We end our brief survey with an example that, once again, involves a continued 
fraction with gamma-distributed elements; the resulting invariant distribution is a 
so-called generalised inverse Gaussian distribution. 
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a; > 0, 



Example 3 (The generalised inverse Gaussian distribution pp. 170-171). Set 

where the a„ are independent, gamma-distributed random variables with parame- 
ters p and s. It was shown by Letac and Seshadri '181 that the probability density 
function of the invariant measure h' is then 

where Kp is the modified Bessel function of order p. As we will show below, the 
Lyapunov exponent can be expressed in terms of modified Bessel functions. The 
details can be found in Section See also ^H] and ^ for some generalisations of 
this example. 

In the papers that form the basis of Examples Q and |3| the authors were con- 
cerned — not with products of random matrices — but rather with the problem of 
determining the distribution of a continued fraction with random coefficients. We 
have already mentioned Dyson's continued fraction 1)1. The continued fraction 
studied by Letac & Seshadri is of the form 

(1-7) X = , 



ai 



1 

02 



03 H 

where, as in Dyson's case, the a„ are independent and gamma-distributed. 

In this paper, we generalise this example to the case where the elements take 
values along a ray in the complex plane. Hence, from now on, unless explicitly 
stated otherwise, we consider matrices in GL(2, C). 

1.1. Main results. Fix a constant 

a e (-7r/2,7r/2) 

and consider the one-parameter family of complex matrices of the form 

A(a):=(; J. 

where the a„ are, again, independent gamma-distributed random variables. 
The corresponding random continued fraction is 

(1.9) Z = ^—^ . 

aie'" H 



1 

aze"- 



age'" H 

The random variable Z takes values in the cone 

(1.10) Sa = {zeC: |argz| < \a\} . 

Our main contribution is an explicit formula for the distribution of Z or, equiv- 
alently, for the /i-invariant measure v associated with the infinite product of the 
random matrices H1.8|l . 
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(c) (d) 

Figure 1 . The probability density function of Z ior p = s = 1 and 
four different values of a: (a) 7r/20, (b) tt/6, (c) tt/S and (d) 9it/20. 
Blue and red correspond to low and high probability respectively. 
For values of a close to (resp. tt/2) the density is locahsed near 
the positive real half- line (resp. the imaginary axis), see Theorem 
Hand 113 



Theorem 1. Let < |a| < 7r/2. Suppose that the a„ are independent, gamma- 
distributed random variables with parameters p,s > 0, and write z = re'^. Then 
the probability density function of Z — equivalently, that of the ^-invariant measure 
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V — is supported on Sa and given by 
sin(2|a|) 



(1.11) Uz)^ 



2Kp{2e'"/s)f r2 sin2(a- 



-0) 



sin(Q; — 6) 
sm{a + 9) 



X exp 



sin(2a) 



1 



r sin(Q; — 6) sin(a + 6 

In particular, the density is a smooth Junction that decays exponentially fast at 
infinity in every direction contained in the cone S(y . 

Plots of the probability density function for p = s = 1 and various values of 
a are shown in Figure ^ 

In Section 1^1 we use the above result to calculate the corresponding Lyapunov 
exponent and express it in terms of the logarithmic derivative of the modified Bessel 
function, namely (see Theorem 0} 



Xp,s{a) = Re 



2 ia\ 



By considering the weak limit in Hl.ll|l as a — > 0, we recover the distribution 
found originally by Letac and Seshadri (see Section [SJ. 

The behaviour of the random variable Z as |a| — > 7r/2— is particularly interest- 
ing. Figure |2] shows plots of the probability density function of arg Z for various 
values of a when p = s = \. As |a| approaches 7r/2, the support of the measure v 
becomes concentrated on the imaginary axis. The weak limit as \a\ — s- 7r/2— leads 




Figure 2. The probability density function of argZ for p = s = 1 
and various values of a. 

to our second result: 

Theorem 2. Let \a\ = 7r/2. Suppose that the a„ are independent and gamma- 
distributed, with parameters p,s > 0, and write z — x + iy. Then, the probability 
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density function of Z- 
given by 

f±7r/2{z) = 



- equivalently, that of the ^-invariant distribution v — is 
S{x) 



[j^{2/s) + Y^^2/s)] 



1 

P+i 



exp 



y 



exp 



y 



whe 



4y) 







ify < 
ify > 



1 



tP-^ dt. 



and S is the Dirac delta. In particular, v is supported on the imaginary axis, where 
the density decays algebraically (like l/j/^J at ±oo. 

Plots of /7r/2(i2/) for various values of p and s are shown in Figure|31 For compar- 
ison, the figure includes a plot of the Cauchy probability density function. We shall 
see in 21 that, for a — ±7r/2 and s small, Z is approximately Cauchy-distributed 
along the imaginary axis. 




Figure 3. The probability density function 7,^/2 for x = and 
various values of the parameters s and p: (a) p = 1; (b) s = 1. 
The dashed curve in (a) corresponds to the Cauchy probability 
density function -^j^-jp. 



Next, we discuss some interesting applications of these results. 

1.2. Random Schrodinger operators. In 1958, P. W. Anderson (2j| postulated 
that the spectra of Schrodinger operators with random potentials should exhibit a 
tendency towards a dense pure point spectrum with bound eigenstates. When the 
eigenstates decay exponentially fast, one speaks of strong or exponential localisa- 
tion. This type of localisation has been proved for one-dimensional systems; cf. [7|, 
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Chap. VIII. Highcr-dimcnsional systems wiU not be considered here; we refer the 
reader to |Zj for a discussion of that case. 

In order to make the connection between this general theory and our own results, 
let us first note that by conjugating the matrices (|f .8(1 for a = -k jl with a rotation 
that maps P(R^) into the imaginary line, we obtain 

K* An{oi)R — a^) ' where R~{^ ^ 

Thus the corresponding random matrices reduce in this limiting case to the transfer 
matrices for the discrete stationary one-dimensional Schrodinger equation 

(1.12) {-^ + V)y^Ey 

with y = {y„}„ez, {l^v)n = Vn+i - 22/„ + y„_i (tight-binding, nearest-neighbour 
laplacian), (V y)n — oLnVn (independent identically distributed random potential) 
and energy level E = 2 (cf. 6 , p. 187). 

It is readily seen that the above equation can be rewritten in the form 

{Vn Vn+l) = {Vn-l Vn) (^J = (-i)" (^0 yi) R* Mn{7r /2)R . 

It is well known TT" that the spectrum on ^^(Z) of the Hamiltonian associated 
with Equation (|1.12() . i.e of the operator H = —A + V is given, with probability 
one, by 

(1.13) a(i/)=a(-A)+supp(M), 

where (t(— A) = [0,4] and /i denotes the probability measure of the potential. 

It has been shown in 8': that, as long as fi is not concentrated on a single point 
and possesses some finite moment, the spectrum ct{H) is pure-point with probability 
one. The eigenfunctions are exponentially localised with a rate determined by the 
value of the Lyapunov exponent (see 0in|)- Thus is view of our results we have 

Proposition 1.1. Let H be the hamiltonian of the Anderson tight-binding model 
in dimension 1 with ^p_s- distributed potential. The exponential localization rate for 
the energy level 

E = 2e a{H) = (0,+oo) 

is given by 

dpKp{2i/s) 



A,,p(7r/2) = Re- 



Kpi2t/s) 



We note that the Lyapunov exponent is a smooth function of the energy level 
E (cf. |7], Proposition VIII. 1.1). Thus, although the above result is not sufficient 
to give estimates of the decay rate of all the eigenstates, it does yield quantitative 
information in the close enough neighbourhood of i? = 2, which is almost surely in 
the bulk of the spectrum. 



1.3. Random Stieltjes functions. Another motivation for this work is to study 
the convergence of Pade approximation for series that are random in some sense 
Consider the class of random Stieltjes functions whose continued fraction 
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expansion is 

1 



Cl 



t 

C2 



C3 + • • • 

where the c„ are independent draws from the gamma distribution with parameters 
p and a. By truncating this continued fraction, we obtain rational approximations 
of i^: 

(1.14) Fn{t):= , teC\R-. 

+ ^ 

C2 H \ 

It is a well-known fact that the Fn are diagonal or near-diagonal Pade approximants 
of F Pl ^ . Roughly speaking, this means that 

F(t) - F„(<) = 0(t") ast->0, 

and that there is no other rational function — with a numerator and a denomina- 
tor of lesser degree — with this property. The rate of decay of the error of Pade 
approximation as n — > cx), t fixed, is of considerable practical interest. 

A simple calculation reveals that the distribution of ^/tF{t) is precisely that of 
the random variable Z defined by Equation H1.9|l if we take 

c„ cr argi 
On ~ — =, s ~ — ;= and a = . 

The rate of convergence of Pade approximation for a typical realisation of the 
function F is then given by 

. _2A, _ almost surely. 

n n^oo '^'v^ \ 2 / 

The proof of this result requires a careful study of the ergodic properties of the 
Markov chain associated with the continued fraction; a detailed treatment will be 
given in a separate publication. 

The remainder of the paper provides the details of the proofs of these assertions. 
In the next section, we introduce some notation, elaborate the correspondence 
between products of 2 x 2 matrices and continued fractions, and derive an integral 
equation for the unknown probability density function of Z . In [21 we deduce a 
partial differential equation, and solve it. The case a = ±7r/2 is treated in SQl^nd, 
in 121 we show that the singular measure obtained in this case is a weak limit of that 
found earlier. Finally, fJHlis concerned with the explicit calculation of the Lyapunov 
exponent. 

2. Preliminaries 

The correspondence between continued fractions and products of 2 x 2 matrices 
is well-known (see for instance 0, Part A, §VI.5), but since it forms the basis of our 
approach, it is desirable to elaborate it here. The starting point is the observation 
that any linear fractional transformation 

aiiz + ai2 



Hz) 



0212 + a22 
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corresponds to the action of a matrix 
(2.1) A = 



an oi2 

on a complex or real projective line P(C^). Indeed, denoting by 

-ft 

any nonzero vector in and defining its projection by 

Uz 



oo otherwise ' 



we obtain a natural identification of P(C^) with IIC^ = C = C U {oo}. The action 
of A on P(C^) (denoted by •) can then be defined via 

(2.2) a-z = ua(';)=^^^^^±^ = h^) 

Then, for the choice 



I J a2iz + a22 



(2.3) A^Aia):=(^^ , |a| < |, 
we get 

(2.4) A-z^—^=:J^a{z). 

z + ae'" 

It is convenient to extend the notation for the set Sa, defined by H1.10|l for 
< |q;| < tt/2, to the limiting cases a = 0, ±tt/2. We shall write 

Sq := M+ and S_^/2 = 5*^/2 := iK = {z G C : | arg z| = 7r/2} . 

Then, for every a > 0, 

Now, suppose that a is a positive random variable distributed according to a prob- 
ability measure /i(da) on M+, and let Zq be an arbitrary random variable taking 
values in Sa- Consider the iteration 



(2.5) Z„ = J^a,, i^Zn-i) , n e N, 

where the a„ are independent draws from (M_|_ , /i) . Under appropriate conditions 
on ij, 12], it can be shown that the Markov chain defined by Equation (|2.5() has a 
stationary distribution i/(dz) independent of the starting value Zq. Moreover, the 
latter distribution coincides with the distribution of the random continued fraction 
(|1.9|) . which can be constructed by means of the successive applications of the 
corresponding random matrices 

Mnia) ■ z Ai{a)A2ia) ■ ■ ■ An{a) ■ z = Ta^ o Tao_ o • • • o TaS^) 



aie" 



a2e'" + 
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Hence, the /x-invariant measure on the projective Hne associated with the infinite 
product of the An is precisely p{dz). 

Let us denote by Vn the distribution of Z„. We are only interested in the case 
where the measures t'„(dz) and u{dz) on Sa are absolutely continuous with respect 
to the Lebesgue measure d^; so let and be the probability density functions 
of Zn and Z respectively, i.e. 

i/{dz) = fa{z)dz and ^'^(dz) = /„^„(2;)d2 for z £ . 

Our aim is to find an explicit formula for /„. To this end, we derive a recurrence 
relation for the and hence, by taking the limit as n — > oo, find an integral 

equation satisfied by f^. We shall assume that „ and fa are smooth. We 
introduce the map Qa '• Sqc — ^ • — ^ defined, by 

Qaiz) = -- ae'« . 
z 

Thus, Qa is the left inverse of Ta- 

yzeSc, {Ga ° :Fa) {Z) = Z . 

On the other hand, for z £ Sa, 

(^a O Qa) (Z) = Z 

only if Gaiz) G Sa- 

Let 5 be a measurable set. We have 

Pr (4 G '5) = y Pr (Z„ G 5 I a„ = a) fi{da) 

= / Pr (j^a{Zn-i) e S) fi{da) = / Pr e GaiS) n Sa) M(da) 

ja,n-\{.z)dz\i{da) 

JAJga{S)nSa 
where 

^:={aeR+ : ea(5)n5„ ^0} . 

For a G A, we have 



L 



fa.n-l{z)dz = / fa,n-l{Qa{w)) Jaiw)dw, 

where jTa is the jacobian of the transformation 

Z = Ga{w) . 

Hence 



Pr Z„e5 =/ / fa,„-l{Ga{w)) Ja{w)dwil{da). 

Now, consider the region 

R := {{a, w) eR+x Sa-. aG A and w e J^a {Ga{S) D Sa)} ■ 
It is readily verified that 

R = {{a,w) eR+ X Sa ■■ w e S and a e A{w)} , 

where 

(2.6) A{z) := {a e IR+ : Ga{z) S Sa} • 
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Hence, 

Pr (Zn eS) ^ [ [ U.n-l (Gaiw)) Ja{w) /i(da) dw . 

^ ' jsJa(w) 



It follows that 

(2.7) /,,„(z) = / {Gaiz)) Ja{z) ^x{Aa) . 

Ja(z) 

Taking the limit as n — )■ oo, we conclude that, if the continued fraction H1.9|l has 
a smooth probability density function /„, then fa satisfies the equation 



(2.8) fa{z)^ faiGaiz)) Jaiz)f^{da) V Z G 5^ . 

JA{z) 

We shall henceforth assume that 

/i(da) = jp,s{a)da . 
The Lyapunov exponent will be denoted by 

(2.9) Ap,,(a) := lim -\n\Mn{a)\ ■ 

n^oc n 

Here, 

\M\ = sup \Mz\ , 

|z| = l 

where | • | denotes the standard euclidean norm on C^. We can express the Lyapunov 
exponent in terms of the measure v as follows 



(2.10) Ap,,(«)--/ ln|z|i/(dz) = - / ln|z|/,(z)dz. 
A proof is given Appendix lUl 

Remark 2.1. An equivalent way of presenting our results is to work with the recip- 
rocal of the continued fraction (|1.9(l . i.e. 

(2.11) aie'" + ^ 



1 

aze"" 



ase'" H 

The corresponding linear fractional transformation is 

z 1-^ ae'" + 1/z 
and so the matrices in the product (|l.l|l are of the form 



(2.12) 



a&" 1 
1 



In particular, when a = tt/2, we have 



ae'" l\ fan — 1 



Readers familiar with the application of the theory of products of random matrices 
to the spectral theory of Schrodinger operators will easily recognize the Schrodinger 
transfer matrix on the right-hand side of this equation. 
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Finally, we remark that, if < a < 7r/2 and fa is the probability density function 
of the random continued fraction H1.9() . then the probability density function, say 
u, of (|2.11|) is given by 

1 

It will sometimes be convenient to work with u instead of fa- Note that in this 
case, i.e. for the matrices of the form H2.12(l the Lyapunov exponent Ap,s takes the 
form (see Appendix 

(2.13) Ap,^(a) = / ln|z|w(z)dz. 

Letac & Seshadri 18 considered the particular case a — 0. By using the Laplace 
transform, they showed that the density function of the stationary distribution is 
given by Equation H1.6(l . In this paper, we study the case < a < 7r/2 by a 
different method. Our approach is inspired by Alain Comtet's derivation of (|1.6|l 
in the particular case a = and p ~ 1 jlOj : By using the fact that 71, ^ is a 
simple exponential, he showed how to replace Equation (|2.8|) by one that involves 
fo{l/x), fo{x) and a derivative of fo- We shall develop this idea and seek to derive 
a differential equation for the unknown density function fa , making use of the fact 
that, for p € N, the probability density function 7^ ^ of the gamma distribution 
solves the differential equation 



— + - 

da s 



" n ^!l2lhl(^^^ JO if < /c < p - 1 



(2.14) 

The key observation is that this is a linear equation with constant coefficients. 



3. The case < |a| < tt/2 
Equation (|2.8|) states that 

(3.1) faiz) ^J^J^ ^" " ^"^^ j ^ ^ 
where 

(3.2) Nsin(a + argz)^ 

sin(2Q;) 

For convenience, we use the notation 

Rez —: X, Imz —: y, \z\ —: r and argz —: 9 

and set 

(3.3) u(a;,y) := -^/a(l/z). 
By Equation H3.1|l . we then have 

pa{x,y) 

(3.4) u{x,y)^ / fa{z-ae'°')-fp^sia)da 

Jo 

where 

X sin a + y cos a 



a 



sin(2a) 
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For the sake of clarity, let us begin by considering the case p = I. In order 
to exploit the fact that 71 _s solves Equation (|2.14(l . we introduce the first-order 
differential operator 

(3.5) 80.:= cos a— — hsina— 

ox ay 

and apply it to both sides of Equation (|3.4|l . Using 

daa{x,y) = 1, 

we obtain 

pa{x.y) 

dau{x,y) = Ja{z - a{x,y)e'°')^i^s{a{x,y)) + / (9a/a(z - ae'")7i,s(a) da 



(•a{x,y) J 

= /a(z - a(a;,y)e'")7i,s(a(x,?/)) - / — /^(z - ae'")7i,s(a) da 



Hence, we have 

(3.6) dau{x,y) + -u{x,y) ^ -^u ' ^ ^ 



Now, the characteristics of the first-order operator da are the curves of equation 

z(i) = z(0) + ie'". 

Thus, the functions that are constant along the characteristics depend only on 
ycosck — a; sin a. This suggests the introduction of the new independent variables 
ro V xsina + ycosa 

(3.7) X= 5 ^ and 1 = ycosa — xsma . 

X + y 

Set 

(3.8) U{X,Y) ■.^uix,y). 
Then, Equation (EU gives 

BU 

Evaluating dU/dX at {Y, X) instead of {X, Y), we deduce that 



§(r,x).^[.M^,y)-L/(r,x)] 



[r-^U{Y,X)-U{X,Y)] 



2xys 

Hence, we have the symmetry property 
(3.10, ^ 

This property can be used to eliminate the term U(Y, X) from Equation (|3.9|l but, 
although the resulting second-order partial differential equation for U is linear, the 
presence of coefficients makes it awkward to solve. We shall seek to guess the form 
of the solution from Equation (|3.9|) instead. 
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The first observation is that, when 1/s = 0, the solution is independent of X. 
Secondly, we know that the solution enjoys the symmetry property 1)^11 Hence, 
we are led to the ansatz 

(3.11) U{X,Y) = C<^'(Y) exp |i + $(F)]| , 

where $ is some univariate function to be determined, and the prime symbol denotes 
differentiation. Equation (|3.9() then gives 



Set 
Then 

or, in terms of x and y 

I 2 I 2\2 

[x +y ) (p 
Let ip{t) = ct^. Then 



r^ip{X) = (^(F) - 2xy 



X sm a + y cos a 
x^ + y"^ 



ip{y cos a — X sin a) — 2xy . 



c [x sin a + 2xy sin a cos a + y cos a) 
and so 



/ 2 2 2 2 \ 

c (x sin a — 2a;?/ sin a cos a + y cos a) — 2xy 



Hence, we have found 



and so 



sin(2a) 



sin(2a) 

, $'(t) = 
sin(2a) 



sin(2Q;) 
i2 



Reporting this in Equation (|3.11|) we obtain 
UiX,Y)=^c'^e.p 



sin(2a) 



1 1 

X^Y 



Hence, in terms of /, r and 0, we have found that, when a is 7p_s-distributed with 
p — I, then 

sin(2a) 



sm [a + 0) 



r sin(Q; — 9) sm{a + 9) 



where Ci^s{ol) > is a normalisation constant. 

For p = 2, 3, . . ., we find by a similar calculation that the equation for u is 
IIP 



(3.12) 



da 



u{x,y) 



1 



-y 



2-2 _j_ j^2 ' ^2 _|_ y2 



It would be tedious to analyse this equation in detail. It is however natural to 
conjecture that, as in the case a — 0, u consists of an exponential term independent 
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of p, multiplied by some algebraic part independent of s. Thus, we seek a solution 
of the form 



(3.13) u{x, y) = gp{x, y) exp | '^^"(^^"^ 

By setting 1/s — in Equation (|3.12|) . we find 

da9p{x^y) = 0. 
The general solution of this equation is 



1 



1 

Y 



gp{x,y)^MY)+r^X^i{Y) + ... + {r^Xy \j>^i{Y), 

for some univariate functions tpk, < k < p. By substituting the resulting expres- 
sion for u in Equation (|3.12|) . we obtain 



Ct-P~^ if/c=p-l 
otherwise 



These results are summarised and stated in the most general form in Theorem 
n The proof is a simple calculation carried out in Appendix IdI 

In Appendix ^ we show how to express the moments of the distribution ly in 
terms of products of modified Bessel functions. From these calculations, we deduce 
in particular that 

(3.14) Cp,s{a) = i ^ . 

|2i^p(2e-/s)|' 

Furthermore, the mean and the variance are given respectively by 

Kp^i (2e-'"/s) 



(3.15) 
and 

(3.16) Var(Z) 



E(Z) 



z iy{dz) 



iE(z)r 



Kp (2c--/s) 



-Im 



^i, Xp_i(2e-'"/g) 



sin(2a) 

4. The case a = ±tt/2 

For \a\ = tt/2, we introduce a function f± : M ^ M defined by 

(4.1) f±^/2{z) = Cp,s{±n/2) f±{y)5{x) for z ^ x + iy with a; > and y e K, 

where Cp,s(±7r/2) is the normalisation constant and S is the Dirac delta. In order 
to find an equation for f±, we write the iteration (|2.5|l in terms of Y„ G M, where 

Zn = ^Yn ■ 

This gives 



(4.2) 



Yn = 



1 



±a„ + Yn-i 



n e 



It is then straightforward to obtain 

(4.3) f±{y) = \rf± 

V Jq 

Since 

(4.4) VyeM, 



y 



fi{da), 
/_(y) = /+(-y), 



ye 
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we need only consider the case a — tt/2. 

As in the previous section, we shaU seek to replace H4.3II by a differential equation. 
Set 

u{y) 



Then, Equation (|4.3() gives 



y 



f+iy - a)7p,s(a) da. 



Differentiating with respect to y, we obtain 



du , 
dy 



iy) = - 



da 



f+iy - a)7p,s(a) da 



-f+iy - ahp,sia)\^ + I f+iy ~ a) (a) da 



= U ivhpAO) + f+iy ~ (a) da . 

By using the fact that 7p^s satisfies (|2.14|) and repeating this calculation, we readily 
deduce the equation 

p 



(4.5) 



uiy) 



— + - 

dy s_ 

Again, in order to solve this equation, it is useful to consider the case p = I first. 



We then have 
(4.6) 

and so we deduce the property 
(4.7) 



dw, , 1 
dy s 



ui-l/y) - yuiy) 



1 du du 

--T-i-^ly) = v:riy>- 

y dy dy 



This property can be used to eliminate the m(— 1/y) term in Equation (|4.t)|l . Mul- 
tiplying the equation by y, differentiating, and then multiplying by — y, we obtain 



y 



dy 



y^^iiy) 

dy 



~i-^/y) + y^ [v'^uiy)] 

ydy dy 



Hence, after simplification and integration, we find 

-y 



2dU 1 , o^ ^ 

d^ = I(^+^)"+^ 



for some constant C; integrating again, we obtain 



u{y) — C exp 



1 / I 



s \y 



t exp 



-(- 

s \t 



In terms of /+, this gives, after some simple manipulations, 



1 




r 




^exp 




/ exp 











dt. 



dt 
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where 

^' \0 ify>0 

This result generahses to arbitrary positive values of p; see Theorem |21 which is 
proved in Appendix IeI 

Remark 4.1. The rigorous determination of the normalisation constant requires 
some care, and we shall in fact make use of a weak limit result which will be 
discussed in the next section. For the moment, let us merely point out that the 
function f± thus defined is integrablc. Indeed, it is readily verified (by use of 
L'Hospital's rule) that the function is continuous at y = 0. By construction, it 
satisfies Equation (|4.3|) . Hence 

/±(y)=0(y"^) as|y|^oo. 

This shows that f± e L^(M). We shall see in due course that 

(4.8) Cp,«(V2)= lim^ ^ ^ 



— /2- 4 \Kp (2e-/s)|' [J^iV^) + Yp^2/s)] 



Plots of Cp^s{'^/'2)f+ for various choices of the parameters s and p are shown in 
Figure |31 We can gain some insight into the behaviour of the random variable Z 
by considering the limit s — > (cf. Appendix IF|| . Set 



and note that 



1 = 



t2 + l 



Then, using integration by parts, 



iv) 



exp 



' s\t * 



s exp 



/ -¥>'(*) exp 

Jc{y) S 



y 



,p+i 



i + r 



exp 



1 . ^ 
S 



d_ 

dt 



dt 



Iterating, we readily obtain 

oo 

(4.9) /+(y)^5]/|")(2/)5" ass 

where 



(4.10) /f(y) = 



1 



and fi-+'\y) 



-1 



dt. 



y'^'ft\y) 



l + y2 yP-i(l + y2)dy 

So by taking into account that (cf. §13.74) 

TT^ [j2(2/s) + rj,2(2/s)] -TTs as s -> 0, 

we see that, for s small, Z is approximately Cauchy-distributed along the imaginary 
axis; Figure |31 (a) provides a graphical illustration of this fact. 
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5. Weak limits of the random variable 



Given the Letac-Seshadri result expressed in Equation 1)1.6(1 . the form of the 
density function for a ^ 0, given in Theorem^ may seem somewhat surprising. In 
this section, we show that the random variable Z has a weak limit as a ^ 0+, and 
that the density function of the limit is indeed given by the formula 1(1. 6|l . In the 
same way, it will be shown that Z has weak limits as \a\ — > tt/2— whose density 
functions are those found in ^ In fact, it was by considering this weak limit that 
we were able to discover the formula for /±7r/2 when p > 1. 

For this purpose, it will be convenient to extend the definition of fa to the whole 
of the complex plane as follows: 

Then, we define a probability measure z/q on sets 5 C C via 



U{z)dz, 0<|a|<7r/2. 



We also define the following singular probability measure on 



MS) 



1 



2Kp{2/s) 



w 



X ^ exp 



SnR4 



da; . 



We shall use the notation — > to indicate a weak limit (limit in distribution). 

Theorem 3. Let X , iY± and Z^ he random variables with probability measures vq, 
v±T^/2 o,nd Va respectively. Then 

w , ^ w 



X and Zr. 



Proof. Let g : C ^ M be an arbitrary bounded and continuous function. It will be 
necessary and sufficient (see |^) to show that 



(5.1) 



g{z)fa{z)Az 



g{z)vt{dz), ^ = 0, ±7r/2 =F • 



Furthermore, since 

f-a{z) = fa{z), 

we shall only consider positive values of a. 

We express the integral on the left-hand side in polar coordinates, and then make 
the substitution 

sin(Q! - 



(5.2) 
Note that 



t 



dt 



sm{a + 6) 
r sin(2Q; 



sin {a + I 



> 0. 



-de. 



t + - 



sm t/ = 



cos ( 



r cos a 
sgn(r - t) 



sin(2Q!) 1 sin(2Q;) 



sin(a + 0) r sin(Q; 







if ^ = 0- 
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and 



1 



cosy = 



So, the substitution H5.2|l yields 



^ I 1 if ^ = 0+ 

o^-^i \ if^ = 7r/2- 



(5.3) 
where 



oo roo 



"'0 



(fa {r, t)dMr, 



tP- 



(Pa-= go z{r, t) cxp 



-1 / cos^ a 



r + t 



sin a [ t + — + r 



and 



Cp,s{a) 



1 



|2ifp(2e-/s)|^ 

Our intention is to let a tend to its limit and take the limit on the right-hand 
side of Equation H5.3|l under the integral sign. For this purpose, we need to find an 
integrable function that provides an upper bound for the family {va}- By using 
the elementary inequality 



we obtain 
where (p : 



r + 1 



is defined by 



sin a > min < 1, 

^ ' r + t 



<Va < ip, 



l\r-t\ / 1 1 
s r + 1 \ t r 



\r-t\\ \r-t\ 



r + t 



l5lU~(C) 



Furthermore, as will become clear very shortly when we consider the limit a 
7r/2— , we have 



I5l 



L~(C) 



oo poo 



ip{r, t)drdt 



exp 



exp 



c(y) 



' s\t * 



Cp^,(7r/2)-i 



We conclude (see Remark l4.1|) that ip is integrable and so we can apply Lebesgue's 
Dominated Convergence Theorem: 



(5.4) lim / g{z)va{dz) = lim Cp^s{c() lim 



oo pCO 



oo /"OO 



^0 



ipa{r, t)dtdr 



poo />oo 

lim Cp s{ce) / / lim ipa{r,t)dtdr 
="^0+ ' Jo Jo 



OO /'OO 



AK^p{2/s)Jo Jo ' rP+^ 



.9(^) 7;;xTexp 



1 1 

hr + - 

t r 



\dtdr 



PRODUCTS OF RANDOM MATRICES 



21 



The double integral on the right-hand side gives 



r P ^ exp 



1 



(r + 1/r) > g(r)dr 



By using the substitution u = 1/t, and reporting the result in Equation H5.4|l . we 
obtain 



(5.5) ^lirn^ / g{z)iyaidz) 



Similarly, 



g{z)iyQ{dz) 



(5.6) lim / g{z)i>a{dz) = lim Cp,s{a) lim 

a — >7t/2 



oo roo 



lim Cp^s{a) 



a — *7r/2— a — >tt/2— Jq Jq 

oo POO 

/ lim ipa {r, t) dtdr 

^0 a^7r/2- 
1 



(PaiT, t) dtdr 



oo /"OO 



TT^ [j2(2/.) + r/(2/s)] 
l|r-i| 



5(sgn(r-t)r)- 



s r + i 



1 1 

t H hr+ - 



didr . 



The double integral on the right-hand side equals 



oo />r 



7-P+l 



exp 



1 1 

r -t^ 

t r 



dtdr 



oo roQ 



9{~r) exp 



7-P+l 
1 



Jr 
oo />'r 

g(r) t^^^ exp 

Jo 

" '"■ /I 

5(r) (-t)''-icxp( - 



t - r + 



1 1 



r t 



\dtdr 



CX2 ^ — OO 



t 

1 

t - - 
t 



9{y) ■ 



1 



exp 



1 

y-- 
y. 



dt- 

dt- 

y 

iv) 



rP+l 



exp 





I 






r 


) 




r 





dr 



exp 



(-r)P+i 
r iP-iexp(i 



1 

s 

1 

t — 



1 

r 

r 



dtdy 



= Cp.s{tt/2) ^ l_g{z)v^/2{dz) 

Hence, after reporting this in Equation (|5.t)|) . we obtain 

Cp,,(V2)-i 



lim 



.g(z)j/a(dz) 



7r2 [j2(2/.) + yp2(2/,)] 7^ 
By setting = 1 in that equation, we deduce that 



giz)v7r/2{dz) . 



Cp,s{n/2) 



7r2 [j2(2/.) + i;2(2/,)] 



and the proof is complete. 



□ 
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6. Calculation of the Lyapunov exponent 

The Lyapunov exponent associated with the random continued fraction 
features prominently in the apphcations considered in the introduction. By using 
the formulae given in Theorems ^ and [2 its computation reduces to a problem of 
quadrature. In this section, we prove that the Lyapunov exponent Ap^s(a) is in 
fact given by the logarithmic derivative of an appropriate modified Bessel function 
for arbitrary p, s > and a € [— 7r/2, 7r/2]. In particular, the result yields explicit 
formulae for Xp^s{a) when p € N, and simple asymptotic expansions for small and 
large s when p > 0. 

Theorem 4. For any p, s > and a [— 7r/2, 7r/2], the Lyapunov exponent Xp^sict) 
associated with the invariant measure v takes the form 

(6.1) XpAa)-Re^p^fi^. 

Proof. Using equation (|2.1U|) and the substitution 
sin(a — 9) r- , , 

. ) , „ , p^rVt, At)=t + 1 t + 2cos{2a), 

we have 

|2ifp(2e'7s)|'Ap,,(a)--|2i^p(2e'7s)|' / ln|z|/„(z)dz 



OO /"OO 



— In t — In p 



^-^exp|-^(p+l/p)l^dt 



s p 



lnt.-^i fexp -^(p+l/p)U^dt. 



2 







Now using the relation 



dp 

together with the uniform convergence (for p in compact sets) of the above integral 
and reversing the change of variables we can continue as follows 



\2Kp{2e^ys)\'XpAa) = \^J^ ^p^i exp |-^(p + l/p)| ^dt 
1 A \2Kp (2e'"/s) f = 4Re i^Kp (2e'"/s) Kp (2e-'"/s) 



2 9p I " ^ ' \dp 

Hence, 

Re (j^Kp (2e-/s) Kp (2e--/. 



|Xp(2e-/s)r 



-^Kp{2e'°'ls) d (2 ■ 

Re ' ^Re^lnKp -c^ 

Kp (2e'"/s) dp ^ 



□ 
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In the remarks below, we mention some of the most important consequences of 
the above theorem. 

Remark 6.1. It is easy to derive a general recurrence relation for the Lyapunov 
exponent. To this end, it is useful to define its complex generalisation as follows: 

(6.2) A,H:=9,lnif,H = ^^. 

Kp{w) 

Differentiating with respect to p the following well-known recurrence relation for 
the modified Bessel function (see [21 §3-71) 

one obtains the recurrence relation for Ap{w) 

which, in view of the fact that 

Ap,,(a) =ReAp 0e'" 
provides an efficient means of computing Ap_s(Q;). 




Figure 4. The Lyapunov exponent Xp_s{ct) as a function of s and 
a for p = 1, 2, 5. 
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Remark 6.2. When p = n e N, we have 

n-l 



(6.4) 



7 I 



^ 2(n- fc)fc! V s 



fc=0 



k—n 



Indeed, the result follows from Theorem 0] and the formula (see §9.6.45 in ^J) 



d_ 

dp 



Kp{w) 



n-l 

E 



2{n-k)k\ 



-I P=" fc=0 

The two most important special cases of (|6.4(l are: 
(1) The Lyapunov exponent for the generalized inverse Gaussian law. 



n-l 



(6.5) 



A„,.(0) = 



k=0 



s"-'=n! Kk{2/s) 
2{n-k)k\Kn{2/sy 



(2) The Lyapunov exponent for the Schrodinger case. Using the formula 
(see m §§3.6, 3.7) 

K^{iw) = |e(3»+i)¥(j„(^i;) +iy„(u;)) 
one immediately obtains 



(6.6) A„,,(7r/2) = Y 



k=0 



s'--^n\ Jfc(2/s) J„(2/s) + n(2/g)y„(2/g) 
2{n-k)k\ J2(2/s) + y,?(2/s) 



Next, we investigate the asymptotic behaviour of \p s{oi) for arbitrary positive 
p. In what follows, ip will denote the digamma function. We start with the large s 
asymptotics. 

Proposition 6.1. For arbitrary p > and a G [— 7r/2,7r/2] we have the following 
large s asymptotics of the Lyapunov exponent 



(6.7) 

where 



(6.8) Rp,sia) = < 



Ap,s(a) = Ins + ij{p) + Rp^s{a), 

'2cos(2pa)fJi^if# + 0(l/s2), z/pe(0,l) 
2 cos(2a) + O (In s/s^) , p = 1 



.TSF^ + 0(i?(.)). 



ifp > 1 



{ In s/s^P i/pe(l,2) 
(lns)Vs4 j/p^2 
z/p > 2 

Remark 6.3. Although only the leading term of the expansion for the remainder 
Rp,s is reported here, it will be clear from the proof that Formula (|6.11(l allows one 
to compute arbitrarily many terms. 



PRODUCTS OF RANDOM MATRICES 



25 



Proof. First, we note that, when p e N, the Lyapunov exponent is given exphcitly 
by Equation H6.4|l , and its asymptotics can be deduced directly from the asymptotic 
expansion of the modified Bessel function of integral order (see §8.446 in 

So, from now on, we shall assume that p is not an integer. As a starting point, 
we take the following representations (Equations (8.485) and (8.445) in 15 ) 

(6.9) = - i-pM-ipM 

^ ' ' 2 sin(j57r) 

(6.10) /,(.) ^ {^YtMv){^)'\ 

k=0 

where 

A, ■ L 

"^P- kWip + k + 1) 

These representations yield (cf. Equations (9.6.42-3) in ^) 

dpKp{w) dpi-piw) - dplp{w) 

r. I ^ ^ -7rcot(p7r) + \% . Y(\ — ' 

Kp[w) I-p[z) - Ip{w) 



fc=0 

oo ^, 

= in(-)/p(u;)-(-) E^^-(p)U, 

fe=0 



where 



_ V^p + fc + i) 
''^ fc!r(p + fc + i) 



Thus 

h.p(w) = — 7rcot(p7r) + 

EZo {-Ak,-r> In (f ) + B.^^p) (f - (A^,, In (f ) - B,,p) (f ) 

For small w, we then have the asymptotic expansion 
(6.11) Ap(u;) -7rcot(p7r) +r(l -p)x 

t [a,.^-, 1„ (1) + (I) " + (a,., 1. (i) + (I) 

1 + f f f r(i - p)^,, (I) - f r(i - (I)"' " 

n=l \fc=0 fe=l y 

Applying the identities (cf. Equation (8.365), Points (1) and (8), |15') 

ip{2-p) - tpil-p) = {l-p)^\ -0(1 -P) - 7rcot(p7r) = ip{p), 

and keeping only the most significant terms we obtain, for all positive noninteger 
p and small w, the expansion 

Apiw) ^ H2/w) + ^{p) + 2^!|[^ H2/w){w/2f'' + ^^^^{w/2f + ■■■ 
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The behavior of the error term depends on the value of p. For noninteger p we 
have the foUowing estimates of the error: 0{w^p) for p £ (0,1); 0((lns)/w^P) for 
p€ (1,2) and 0{w^) for p > 2. 

Substituting w = 2e'"/s and taking the real part (of. Theorem^)) completes the 
proof. □ 

It is worth noting that, when p is large, there is an alternative way of obtaining 
the asymptotics of the Lyapunov exponent. Indeed, by using the integral represen- 
tation 

one easily deduces the following version of Equation (8.446) in jl5| : 

^ ' 0<k<p ^ ' 

This yields H6.7|l for p > 2 and explains the absence of "logarithmic" prefactors in 
the low order terms of the expansion for the Lyapunov exponent when p is large. 
Next, we turn to the case of small s. 

Proposition 6.2. For arbitrary p > and a £ [— 7r/2, 7r/2] , we have the following 
small s asymptotics of the Lyapunov exponent: 

oo 

(6.13) Ap,,(a) - Z„s" as s^O, 
where 

, , pcosa , ncos2Q; , »(4»^ — 13) COS 3a 

(6.14) h^- , h^-- , h = -^-^ , 

\ J 1 ^ , 2 8 ' ^ 192 

p(4p2 _ 7) cos 4a , p{-920p^ + A8p* + 1187) cos 5a 



128 ' " 20480 

Remark 6.4. Higher order terms can be derived using formula (|6.16|) below. 

Remark 6.5. When a = ±Tr/2 (and only in this case) the real part of all odd 
order terms in (|6.13|) vanish. Thus the leading order term in this particular case is 
quadratic (cf. Figures ^ and . 

Proof. For small s, the asymptotics follow from the behaviour of the modified Bessel 
functions for large w. We have (see §7.23) 

where 

r(p + TO+l/2) _ (4p2 _ l)(4p2 _ 32)_ ^(4p2 _ (2to- 1)2) 

^' '~ mlT{p-m+ 1/2) ^ 2^"'ml 
and where the constant 8 satisfies 

|e| < 1 and Ree > for p e M and Rcw > . 
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(lfn>p- 1/2, then < 9 < 1, see |211 §7.30). Thus, we have 

^ „ . ^ v^n-l dp{p,m} dp(0(p){p,n)) 
OpKp[w) _ l^m=0 (2u;)"' ^ {2w)^ 



(6.15) 
where 



dp{p,m) 



r(p + m+ 1/2) 



n— 1 (p,m) 
(2^1)^ 



e(p) 



(2u;)" 



m\T(j) — TO + 1/2) 
This leads to the foUowing expansion for large w 

dp{p,m) dp{Q{p){ja,n)) 



[^^{p + TO + 1/2) -■4j{p + m+ 1/2)]. 



(6.16) Kp{w) 



E 



{2w) 



1 + E -EM-e(^):^''"^ 



N' 



N=l 



- (2z«)^ 



{2wy 



From this expansion, every coefficient in H6.13|l can be computed explicitly. In 
particular, by using the basic properties of the gamma and digamma functions, we 
obtain 



(6.17) A,(u;)^fi-£, 



P___P__ Pj'^p'^ - 13) ^ 
w 2w^ 24^3 

p(V - 7) p(-920p2 



1187) 



as w ^ cxD . 



640^5 

By setting w = 2e'"/s and taking the real part, one obtains Equation H6.14I) . □ 
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Remark 6.6. Although the series H6.13|l is divergent, one can compute its sum by 
using standard summation techniques 0]. The plots of Ap,s(7r/2) against s for 
various half-integer values of p shown in Figure |31 (b) were obtained by using a 
Fade approximant — i.e. a rational function of s that matches the series to 0{s^^). 

Acknowledgments: We thank Brian Winn for a very helpful discussion leading 
to the discovery of Formula H6.1|l . 

Appendix A. Some identities involving products of Bessel functions 

We begin with two integral representations of the modified Bessel function: For 
Rez > 0, we have (cf. 24 , §6.22) 

(A.l) Kp{z) t-P-^ exp {-^(t + 1/t)} dt 

When calculating the moments of we will use integrals of the form 
,A,2, X<".(.,.„:4fo.p{-I-^!±^}Av(")."d. 

where n € Z, and the complex numbers u and v are such that 

|argu|<7r, |arg?;|<7r and | arg(u + w)| < 7r/4 . 
For n — —1, the value is given by Macdonald's formula ,24,, §13.71: 
(A.3) I^-^\u,v)^Kpiu)Kj,{v). 

Differentiate this identity with respect to u, and then multiply it by u to obtain 




after integration by parts. If we interchange u and u, this becomes 

Hence, by adding these two identities, we find 
(A.4) lj°Hu,v) = uKp+i{u)Kp{v) + vKp+i{v)Kp{u) - 2pKp{u)Kp{v) . 
By similar calculations, we obtain easily the recurrence relations 

and 



Kp{^) dr. 



(A.5) 



jin) 



d d 

2n ~ u- t; — 

ou ov 



(A.6) 



{v^ ~ zi2)x("-i) = 



d 



d 



du dv 
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The first of these enables tlie evaluation of the integral (|A.2|I in terms of products 
of modified Bessel functions when n = 0, 1, . . ., whereas the second is appropriate 
when n — —2, —3, . . . (provided that ^ w^). 



Appendix B. Calculation of the moments 
In this section, we show how to calculate the moments 



(B.l) 



of the distribution i/. For simplicity, we shall assume that a > 0. 
We have 



M(™-")(a) = 



The substitution 



a Jo 



X exp 



sin(2Q:) 



sin(a 



' ^m+n^i{m-n)e Cp,^(a) sin(2a) 

-2sin^(a + 6') [sin(a + 6') 



1 



r sin(a — 6) sin(a + I 



rdrdO . 



yields 



and 



where 



Hence 



sin(Q; + 9) 



sin (a + 1 



2sin(2a) 



sy^sm{a - 6) s\a{a + 6) s 
(p{t) = i + l/i + 2cos(2Q;). 



h Jo 



oo poo 



m+n — l^p— 1 



-I m—n 



X exp ■ 



r V t 



ArAt 



Next, we make the change of variable 



2 , ^{t) 



ST \ t 
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to obtain 



(B.2) M^™'")(a)=cp,,(a) 



-j-'p—m— 1 

pm+n+1 " ^ ' V / 



X exp 



drdt 



Cp,s{a) 



rn+n „od 

Jo 
X exp 



4-p—m— 1 

— [t+l/t + 2 cos(2a)]" (e- + ie"'")' 



1 



U+- exp -- 



4cos(2a) 



By using the fact that 

1 

if,(z) = - t«-iexp{-^(i + lA)}dt, 

we can, provided that n < m, express this double integral in terms of the integrals 
(|A.2|I . For instance, when m = n = Q, Equation (|B.2|) yields 



exp 



T 4cos(2a) 
= 2cp,s(a) 



T 4cos(2a) 



exp 



t + 



dt- 



dr 

T 

dr 

T 



This can be expressed as 

1 = Af^"i(a) = 4cp,.(a)4-i) , u = w = 2e'"/s . 

So we deduce the value of the normalisation constant from Equation (|A.3p . In the 
same way, by taking m = 1 and n = 0, we obtain the mean: 

E(Z) = M^lfHa) = 4c,,.(«) {vll,-_l^ {u,v) + ul^^'^) m] . 

Its value is given by Equation H3.15II . For the variance, we use 

Var(Z) = M(;/)(a)-|M(i,'°)(a)|' . 

Equation (|B.2p gives 

and so we deduce the value given in Equation (|3.16() . 

Appendix C. The Formula for the Lyapunov exponent 

In this appendix, we consider matrices distributed according to a measure /i on 
SL{2,C) of the form 

a b 



A 



c d 



with exactly one row being random. In other words, the measure /i is such that 
one or the other of the following conditions holds: 



(I) 



a and b are fixed, c and d are random 
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(II) c and d are fixed, a and b are random . 

Let v be ^-invariant measure on the projective space P(C^) = C. In this case, 
we have (cf. Equation (1) in 0, p. 9) 



(p{A ■ z)dfi{A)di^{z) = / ip{z)dh'{z) 
CJSL{2,C) Jp{C^) 

for every bounded Borel function ip. 

First, let us consider the case where Condition holds. Assume that v has 
some positive moment. Suppose also that /i satisfies the conditions required for the 
existence and positivity of the Lyapunov exponent A (see Theorem 3.6 in [H], p. 
27 — which holds also in the complex domain) . Then 

A = / / ln^lilHdM^)dK^)- 

JcJSL(2.C) 

2JcJSL(2X) \z\'^ + l 

Iff In (1 + 1^ . zr) + In d^(^) dKz) = 

i / ln(l + |zr2) diy{z)+ f ln|a + 5z"^| dv{z) 

f ln(l + |z|-2) Au{z)= f ln|a + 6z-i| dt/(z). 

We note that the splitting of the integrals in this calculation is permissible because 
the assumed existence of a positive moment of fi guarantees the integrability of 
ln(l + |z|-2) and ln(a + b\z\-^). 

The case where Condition ((11}) holds involves a similar calculation; hence 



(C.l) 



'/j. In \a + bz ^\diy{z) if Condition (U) holds 
Jj, In |cz + d\ di'{z) if Condition (|lTll holds 



Now, we apply this result to the matrices H2.3|l considered in Section [3 In this 
case. Condition holds; there is only one random entry, and since it is gamma- 
distributed, the support of fi is not contained in any compact subgroup of GL(2, C). 
Furthermore 

A*„{a) = Ai(-a) 

and hence the invariant measure for yl*(a) is continuous for all a G [— 7r/2, 7r/2]. 
We may therefore use Proposition 3.3 in p. 26, to deduce the uniqueness of the 
invariant measure z/ found in the paper. Since, for all a G [— 7r/2, 7r/2], this measure 
possesses positive moments of all orders less then one, we can apply Formula l|C.l|) 
to obtain H2.1U|I . 

For matrices of the form (|2.12|) and for Schrodinger matrices, it is Condition lUIj) 
that holds, and the formula IjClj) leads to H2.13|l . 
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Appendix D. Proof of Theorem ^ 
Proof. By direct substitution in Equation (|3.4|1 . To facilitate the calculation, define 



Va = \/ (x — acosQ!)^ + (2/ — asinQ;)^, 

(a; — a cos a) sin a + (y — a sin a) cos a 



Xa 



and 



We note that 



We can write 
(D.f) /„(z - ae'")7p,,(a) = CY^-^ exp 



Ya = {y — a sin a) cos a — {x — a cos a) sin a . 
ro = r, Xo = X and 1; = !^. 

sin(2a) 1 



s Y 
sPT{p) 



[rlXa] *"+'^exp{i?(a)/4, 



where 



sin(2a) 
£:(a) = — a. 



A lengthy but straightforward calculation reveals that 

2 

E'(a) 



Y 



riXa 

Hence, by integrating Equation (jPTjl with respect to a, we obtain 



(D.2) 

where 

Noting that 



a{x,y) 



fa{z - ae'")7p,s(a) da = C 



~Y^ 



■ cxp 



sin(2a) 1 



a(x,y) 



spr(p) 



£;'(a)exp{i;(a)/s} da. 



_d 
da 



rlXa 



a; sin a + y cos a r^X r^X , 



r^X^ 



4^2 



Y2 



E'ia), 



we obtain, after integration by parts. 



r^X 



p-1 



exp {E{a)/s} 



a{x,y) 



ifp=l 

/exp{£;(0)/s} ifp=l 
1* 



Thus, 



f — 



r^2;^lP-l 



y2 



exp 



sin(2Q:) 1 
s A 
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and, reporting this in Equation (|D.2p . we have shown that 



fa{z - ae'")7p,s(a) da 



C /r^X^"-' 



It follows easily that satisfies Equation H3.4|l 



exp 



sin(2a) 



1 1 



□ 



Appendix E. Proof of Theorem El 

Proof. We need only consider the case a — tt/2. As in the case \a\ < 7r/2, the proof 
is by direct substitution of the expression 



f+iy) 



1 



i/p+i 



cxp 



s \y 



exp 



civ) 



^ s\t * 



tp~^ dt , 



where 



ciy) = 



—oo if y < 
if 2/ > 



into the integral equation (|4.3|l . 

By making the substitution u — —1/y — a in the integral on the right-hand side 
of Equation 14.3|l . we obtain 



(E.l) 











/o \ y 


















— oo 



exp 



1/11 1 
-[---+t+- 

s \u t y 



1 



1 



spr{p) V y 



p-i 



dMi 



At this point, it is helpful to treat the cases y < and y > separately. 

Let us begin by assuming that y > 0. Then c{u) — —oo on the right-hand side 
of the last equation and, by changing the order of integration, we obtain 



1 

Jq 



1 



n.p+l 



exp 



s \ t 



t spr{p) \ u 



p-i 
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Then, after making the substitution t = — followed hy v = 1/u + y, this 
becomes 



1 











Jo 





a ) 7p,s(a)da 



exp 







rv 


1 






j exp 










s 




Jo 


sPT{p) 


-le^ 



where 
(E.2) 



In order to prove that Equation H4.8|l holds for y > 0, we therefore only need to 
verify that 



exp 



s\t * 



tP-^dt= / exp 
Jo 



1 



tP- 



holds identically for every y > 0. Clearly, it holds when y = 0+, and so it will be 
sufficient to show that the derivatives are the same. The derivative of the right-hand 
side is 



exp 



f /I 



s \y 

rv 



- y 

exp 



1 



•iP+ 



1 / I 



exp 



1 /I 



s \ t 



1 



tP- 



jv'{y~t)dt 



tP+i sPT{p) 



(y-tf-iexp 



-iy-t) 



dt 



1 



sPT{p) \t 



1 1 



p-i 



— y^ ^ exp 



1 



5Pr(p) 



iP-^e-'^du, 



after making the substitution u = 1/t — 1/y. The result follows. 

To complete the proof, we also need to consider the case y < 0. Returning to 
Equation (|E.ip . we need to split the range of integration of the variable u into pos- 
itive and negative values. Recalling the definition of c{u) and, as before, changing 
the order of integration, we obtain 
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After making the substitution v = —1/u — y, this becomes 



exp 



s (2/ ^. 



-P+i 



exp 



-i/y 

exp 

-00 



X 



exp 



' s\t * 



where 

^!!>(u;) = 

Therefore, we have to prove that 

rv 



sPT{p) 



yP ^e " dw . 



/' 

J —I 



exp 



' s\t * 



exp 



' s\t * 



I 



exp 



s\t * 



holds for every y < . On both sides, the Umit as y — > —00 is zero. Furthermore, 
the derivative of the right-hand side is 



exp 



1 /I 



s \y 



y 



(_l)p-i 

^—Jtt- im - 1] 



exp 



s [t * 



tP-^cP'i-l/t-y) 



-i/j/ 



exp 



' s\t * 



1 



1 



sPT{p) \ t 
1/!/ 



p-i 



J-00 sPTip) V y . 



- - y ) exp 

p-i 



di 



e = dt 



= yP ^ exp 



1 /I 



s \y 



-y 



1 



s^r(p) 



after making the substitution v = —1/y — t. 

Appendix F. An asymptotic expansion for /+ as s ^ 
It is actually simpler to work with 



□ 



So let write 



where 



u{y) = ^U{-l/y). 



n=l 

«^"^(y) = 4/+^(-V2/)- 
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Equation (|4.1Q(I then gives 
1 



(F.l) 



and 



,(«+i)(y) = ^ 



,p-l 



A direct calculation of the first few terms reveals that the u'"'' have a very simple 

structure: 

(F.2) 



,(2 



4n-3 (2ri-l) 



and 



(2n) 



=n+l ^ ^ ^ 



By substituting these expressions in Equation (|F.ip . one obtains the recurrence 
relations 



u. 



(2n) 
(2n) 
(2n) 

""471 -1 

(2n+l) 

(2n+l) _ 
'n+2 — 

(2n+l) _ 



(p-l + 2n)42"-i), 

(p-3 + 2fc)4iV^-2(fc-2)4l"-^\ 

-2(4n-3)4'„'r3'\ 

(p + 2n)z.|f;J, 



2 < fc < 4n - 2, 



(p + 2n + 2)u^^^l - (p + 2 + 4n)«(f;J , 
-(p~2 + 2/c)4'") 
+ (6 - p - 4fc)u^i"^ + 2(/j - 2)u''^"l , n + 3</c<4n-l, 
^'n"-"'^ = (6 - P - 16n)4?^ 1 + 2(4n - 2)u^^:%, 



u 



(2n+l) 
4n+l 



with the starting value u'"i^ — 1. We then have 



4n-3 



(F.3) 
and 



/ f+iy)dy^ / «(y)d2;^2Vs2-i V 



(2ri-l) 



4n-3 



k—n 



^ ^ WSJ Tit tj 

(F.4) - / In \y\ /./2(2/)dy = / In |y| u{y)dy ^ 2 ^ s'^' E "^^'^/^fc 

■/K JR , 

as s ^ 0, where 

and Pk = 



oik = 



dy 



The Qffc and Pk satisfy the recurrence relations 



Iny dy 

(TTTF 



ttfe+i " I -'^ ^ 2fc ' ^''^^ ^ '^'^ 



ak + f3k TT 
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